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Abstract. We study the closure problem for continuum balance equations that model mesoscale dynamics 
of large ODE systems. The underlying microscale model consists of classical Newton equations of particle 
dynamics. As a mesoscale model we use the balance equations for spatial averages obtained earlier by a 
number of authors: Murdoch and Bedeaux, Hardy, Noll and others. The momentum balance equation con- 
tains a flux (stress), which is given by an exact function of particle positions and velocities. We propose a 
method for approximating this function by a sequence of operators applied to average density and momen- 
tum. The resulting approximate mesoscopic models are systems in closed form. The closed from property 
allows one to work directly with the mesoscale equaitons without the need to calculate underlying particle 
trajectories, which is useful for modeling and simulation of large particle systems. The proposed closure 
method utilizes the theory of ill-posed problems, in particular iterative regularization methods for solving 
first order linear integral equations. The closed from approximations are obtained in two steps. First, we 
use Landweber regularization to (approximately) reconstruct the intcrpolants of relevant microscale quanti- 
ties from the average density and momentum. Second, these reconstructions are substituted into the exact 
formulas for stress. The developed general theory is then applied to non-linear oscillator chains. We conduct 
a detailed study of the simplest zero-order approximation, and show numerically that it works well as long 
as fluctuations of velocity are nearly constant. 



Key Words: FPU chain, particle chain, oscillator chain, upscaling, model reduction, dimension reduction, 
closure problem, 

1. Introduction 

In a series of papers, p], [2], [3], [4], Murdoch and Bedeaux studied continuum mechanical balance 
equations for mesoscopic space time averages of discrete systems. Earlier work of Irving and Kirkwood 14j, 
Noll [TO], and Hardy [T3] on closely related topics should be also mentioned here. The fluxes in balance 
equations (e. g. stress) are given by exact formulas as functions of particle positions and velocities. This is 
useful for linking microscale dynamics with mesoscale phenomena. However, using these formulas requires a 
complete knowledge of underlying particle dynamics. Since many particle systems of interest have enormous 
size, direct simulation of particle trajectories may be intractable. Consequently, it makes sense to look 
for closed form approximations of fluxes in terms of other mesoscale quantities (e.g., average density and 
velocity), rather than microscopic variables. 

In this paper we address the above closure problem for spatially averaged mesoscale dynamics of large 
size classical particle chains. The design of the method was influenced by the following considerations. 

(1) The quantities of interest are space-time continuum averages, such as density, linear momentum, 
stress, energy and others. This choice of averages is natural because these quantities are experimen- 
tally measurable, and also because of their importance in coupled multiscale simulations involving 
both continuum and discrete models. In addition, by working directly with space-time averages in- 
stead of ensemble averages one can bypass a difficult problem of relating probabilistic and space-time 
averages. 

(2) It is desirable to be able to predict behavior of averages on arbitrary time intervals, no matter 
how short. This perspective comes from PDE problems, where observation time is often arbitrary 
and long time behavior is not of interest. When one tracks an ODE systems on an arbitrary time 
interval, transients may be all that is observed. Therefore, we do not use qualitative theory of ODEs, 
primarily concerned with describing long time features of dynamics. This significantly decreases the 
range of available tools. However, the closure problem for mesoscopic PDEs turned out to be a 
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question that can be still answered in a satisfactory way. The methods developed in this fashion can 
be helpful in situations where long time features are not of interest: modeling transient and short- 
lived phenomena, working with metastable systems, and dealing with problems for which relaxation 
times can be hard to estimate. 

(3) We consider particle systems with initial conditions that either known precisely, or ar least such that 
the possible initial positions and velocities are strongly restricted by available a priori information. 
This is in contrast to statistical mechanics, where uncertainty of initial conditions is a major problem. 
In this regard we note that our approach makes sense for discrete models of solid-fluid continuum 
systems, where the smallest relevant length scale is still much larger than a typical intermolecular 
distance. For other particle systems, our method can be used to run deterministic simulations 
repeatedly, in order to accumulate statistical information about the underlying probability density. 

(4) Because of widespread use of computers in physical and engineering sciences, it is useful to develop 
theories tailored for computation, rather than "paper and pencil" modeling. As far as the closure 
problem is concerned, traditional phenomenological approach to formulating constitutive equations 
can be subsumed by a more general problem of finding a computational closure method. In particular, 
a closure method can be realized as an iterative procedure where one inputs the values of the primary 
variables (e.g. density and velocity) computed at the previous moment of time, and the algorithm 
generates the flux (e.g. stress) at the next moment. Then primary variables are updated using 
mesoscopic balance equations, and the process is repeated. In addition, focusing on computing 
one can obtain unconventional but useful continuum mechanical models. By replacing a simple, 
but possibly crude, Taylor series truncation with an algorithm we make it harder to obtain exact 
solutions. Since such solutions are rarely available even for simple classical systems, (e.g. Navier- 
Stokcs equations), this is not a serious drawback. On the positive side, computational closure 
generally contains an explicit (explicitly computable) link between micro- and mesoscale properties. 

(5) An important potential application of closure is development of fast numerical methods for simulating 
meso-scopic dynamics of particle systems. Mesoscale solvers usually employ coarse meshes with mesh 
size much larger than a typical interparticle distance. Then the averages would be usually given by 
their coarse mesh values, while interpolants of microscale quantities are discretized on a fine scale 
mesh. Consequently, a closure method might consist of two generic blocks: (i) reconstruction on 
mesoscale mesh thereby a coarse approximations of fine scale quantities are obtained from averages; 
and (ii) interpolation of the obtained coarse scale discretizations to fine scale. 

The closure algorithm developed in the paper is based on iterative regularization methods for solving 
first kind integral equations. We observe that primary mesoscale averages are related to the interpolants of 
microscale variables via a linear convolution operator. The kernel of this operator is the "window function" 
used in [T] to generate averages. Such integral operators are usually compact. A compact operator may be 
invertible, but the inverse operator is not continuous. Therefore, the problem of reconstructing microscale 
quantities from given averages is ill-posed. Such problems are well studied in the literature [T2 l \6\ \9\ IT5 l fTH] . 
A particular method used in the paper for inverting convolutions is Landweber iteration jS] , [10] . It is known 
that if the error in the data tends to zero, the Landweber method produces successive approximations 
converging to the exact solution. For the merely bounded data error, convergence is replaced by a stopping 
criterion. This criterion provides the optimal number of iterations needed to approximate the solution with 
the accuracy proportional to the error in the data. As a consequence, our method has desirable feature: one 
can improve the approximation quality at the price of increasing the algorithm complexity. This means that 
predictive capability of the method can be regulated depending on available computing power. 

The paper is organized as follows. In Section 2 we describe a general multi-dimensional microscopic 
model. The equations of motion are classical Newton equations. We limit ourselves to the case of short 
range interaction forces that may be either conservative or dissipative. The scaling of particle masses and 
forces reflects a continuum mechanical perspective, that is a family of particle systems of increasing size 
should represent a hypothetical continuum material. As N —> oo, the total mass of the system should 
remain fixed, and the total particle energy should be either fixed, or at least bounded independent of N. 
Next, we recall the main points of averaging theory of Murdoch-Bedeaux and provide mesoscopic balance 
equations and exact formulas for the stress from In Section 3 we develop integral approximations of 
averages, and describe the use of Landweber iterative regularization for approximate reconstruction. Section 
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4 contains the formulation of the scaled ODE equations of the so-called Fermi-Pasta-Ulam (FPU) chains. 
In Section 5 we derive closed form mesoscopic continuum equations of chain dynamics. The complexity of 
these continuum models increases with the order n of the iterative deconvolution approximation. Section 6 
is devoted to the detailed study of the simplest closed model with n = 0, which we call zero-order closure. 
Essentially, zero-order closure means that the microscopic quantities are replaced by their averages. Such 
an approximation can work well only for systems with small fluctuations. To quantify fluctuation size we 
introduce upscaling temperature and the related notion of quasi-isothermal dynamics. For such dynamics, we 
show how to interpolate averages given by mesoscopic mesh values, in order to initialize approximate particle 
positions and velocities. The interpolation procedure is problem-specific: it conserves microscopic energy 
and preserves quasi-isothermal nature of the dynamics. Section 7 contains the results of computational tests. 
Here we apply our zero-closure algorithm to a Hamiltonian chain with the finite range repulsive potential U, 
decreasing as a power of distance. The results show good agreement of zero-order approximations with the 
exact stress produced by direct simulations with 10000-80000 particles, provided the initial conditions have 
small fluctuations. In our example, the initial conditions are such that the upscaling temperature is nearly 
zero during the observation time. We also demonstrate that increasing fluctuations of initial velocities leads 
to a considerable increase in the approximation error, indicating that higher order closure algorithms should 
be used instead of zero-order closure. Applicability of the zero-order closure is further discussed in Section 
8. Finally, conclusions are provided in Section 9. 

2. MlCROSCALE EQUATIONS AND MESOSCALE SPATIAL AVERAGES 

2.1. Scaled ODE problems. The starting point is the microscale ODE problem. In this paper we shall 
work with classical Newton equations of point particle dynamics. The same equations may arise as dis- 
cretization of the momentum balance equation for continuum systems. Consider a system containing N 3> 1 
identical particles, denoted by Pj. The mass of each particle is where M is the total mass of the system. 
Suppose that during the observation time T, p remain inside a bounded domain fl in E d , where d is the 
physical space dimension, usually 1,2 or 3. The positions qi(t) and velocities Vi(t) of particles satisfy a 
system of ODEs 

(2-1) q % = Vi, 

(2-2) ^Vi = ft + ft*, 

subject to the initial conditions 

(2.3) qM = x h Vi(0)=Vi- 

Here f'f^ denotes external forces, such as gravity and confining forces. The interparticle forces f i — J2j f iji 
where are pair interaction forces which depend on the relative positions and velocities of the respective 
particles. 

We are interested in investigating asymptotic behavior of the system as N — > oo. Thus it is convenient 
to introduce a small parameter 

(2.4) e = N- 1 ^, 

characterzing a typical distance between neighboring particles. As e approaches zero, the number of particles 
goes to infinity, and the distances between neighbors shrink. Consequently, the forces in (2.2 1 should be 
properly scaled. The guiding principle for scaling is to make the energy of the system bounded independent 
of N, as N — > oo. In addition, the energy of the initial conditions should be bounded uniformly in N. 

As an example of scaling, consider forces generated by a finite range potential U and assume that each 
particle interacts with no more than a fixed number of neighbors (this is the case, e.g., for particle chains 
with nearest neighbor interaction, where a particle always interacts with two neighbors). The fixed number 
of interacting neighbors implies that there are about N interacting pairs. Assuming also that the system 
is sufficiently dense, and variations of particle concentrations are not large, we can suppose that a typical 
distance between interacting particles is on the order N~ x / d L = eL. The resulting scaling 

(2-5) fii=-^ ufqi q ' 



makes the potential energy of an isolated system bounded independent of N. Kinetic energy will be under 
control provided the total energy of the initial conditions is bounded independent of N. If exterior forces 
are present, they should be scaled as well. 



Remark. Superficially, the system (2.1|, ( 2.2 1 looks similar to the parameter-dependent ODE systems studied 
in numerous works on ODE time homogenization (see e g. |17j and references therein). In the problem 
under study, e depends on the system dimension N, while in the works on time-homogenization and ODE 
perturbation theory, the system size is usually fixed as e — > 0. 

2.2. Length scales. We introduce the following length scales: 

- macroscopic length scale L = diam(f2); 

- microscopic length scale sL; 

- mesoscopic length scale 77L, 

where < r/ < 1 is a parameter that characterizes spatial mesoscale resolution. This parameter is chosen 
based on the desired accuracy, the computational cost requirements, available information about initial 
conditions and behavior of ODE trajectories etc. 

The computational domain £1 is subdivided into mesoscopic cells Cp, f3 = 1, 2, . . . , B, with the side length 
on the order of r\L. The centers xp of Cp are the nodes of the meso-mesh. The number of unknowns in the 
mesoscopic system will be on the order of B. For computational efficiency, one should have B N . This 
does not mean that r\ is close to one. In fact, it makes sense to keep rj as small as possible in order to have 
an additional asymptotic control over the system behavior. Decreasing rj will in general make computations 
more expensive. 

2.3. Averages and their evolution. To define averages we first select a fast decreasing window function 
ip satisfying J tp(x)dx = 1. There are many possible choices of the window function. In the paper we assume, 
unless otherwise indicated, that ip is a compactly supported, differentiable on the interior of its support, and 
non-negative. Next, define 



ip v (x) =77 d tp 



, 7 L 

Once the window function is chosen, we can evaluate the averages of various continuum mechanical 
variables, following pQ, pQ. The mesoscopic average density and momentum are given by 

M N 

(2-6) ^{t,x) = —^, l (x-q i (t)), 

(2.7) p^(t, x) = ^ - *(*))■ 

The meaning of the above definitions becomes clear if one considers ip = {cd) x( x ); where x is a charac- 
teristic function of the unit ball in M. d , and Cd is the volume of the unit ball. Then 

1 (x-q&) 



9 1 c d r, d N ^ X \ ;/ 

The sum in the right hand side gives the number of particles located within distance rj of x at time t. 
Multiplying by M/N we get the total mass of these particles, and dividing by Cd,r] d (the volume of ry-ball) 
gives the usual particle density. 



Differentiating ( |2.6[ ), (2.7! in t, and using the ODEs (2.11, (2.2) one can obtain Q] exact mesoscopic 
balance equations for all primary variables. For example, for an isolated system with (f'f^ = 0), mass 
conservation and momentum balance equations take the form: 

(2.8) d{f + div(p'V) = 0, 

(2.9) d t {jf'v^) + div (p"tf <g) v*) - divT" = 0. 
The stress T v = Tj c) + Tj int) [4j, where 

(2.10) Tj c) (t, x) = -J2 m(Vi - ^(t, x, )) ® ( Vi - v^x, t))^{x - q t ) 



is the convective stress, and 



(2.11) T^(t,x) {lnt) = fa ® (q 3 - Q t ) [ V> (s(x - q 5 ) + (1 - s)(x - g<)) ds 



is the interaction stress. The summation in (2.111 is over all pairs of particles that interact with each 
other. 

Discretizing balance equations on the mesoscopic mesh yields a discrete system of equations, called the 
meso-system, written for mesh values of pj, {jpv n )fj and Tj. The dimension of the meso-system is much 
smaller than the dimension of the original ODE problem. However, at this stage we still have no compu- 



tational savings, since the meso-system is not closed. This means that mesoscopic fluxes such as (2.101, 



(2.11) are expressed as functions of the microscopic positions and velocities. To find these positions and 



velocities, one has to solve the original microscale system (2.1), (2.2). To achieve computational savings we 
need to replace exact fluxes with approximations that involve only mesoscale quantities. We refer to the 
procedure of generating such approximations as a closure method. This closure-based approach has much 
in common with continuum mechanics. The important difference is that the focus is on computing, rather 
than continuum mechanical style modeling of constitutive equations. 

3. Closure via regularized deconvolutions 

3.1. Outline. Our approach is based on a simple idea: the integral approximations of primary averages 
(such as density and velocity) are related to the corresponding microscopic quantities via convolution with 
the kernel Vv Therefore, given primary variables we can (approximately) recover the microscopic positions 
and velocities by numerically inverting convolution operators. The results are inserted into equations for 
secondary averages (or fluxes), such as stress in the momentum balance. This yields closed form balance 
equations that can be simulated efficiently on the mesoscopic mesh. 

3.2. Integral approximation of discrete averages. To exploit the special structure of primary averages, 
it is convenient to approximate sums such as 

1 N 

(3-1) 9 V = ■^229( v j,9j)' l l } v( x ~ 9j) 



by integrals. Since particle positions q^ are not periodically spaced, (3.1) is not in general a Riemann sum 
for gipri(x — •). To interpret the sum correctly, we introduce interpolants q(t,X),v(t,q) of positions and 



velocities, associated with the microscopic ODE system (2.1), (2.2). At t — these interpolants satisfy 

q(Q,X j ) = q° j , v(0,q(0,Xj)) = v° p 
where Xj, j = 1, 2, . . . , N are points of e-periodic rectangular lattice in Q. At other times, 

q(t,Xj) = «!,(*), v{t,q{t,X 3 )) = Vj (t). 



Then we can rewrite (3.1 1 as 

JV 



(3-2) & = t^t ]T ^5 (v (t, q(t, X,)) , q(t, X ,) - q(t, X,)), 

3=1 



where |Q| denotes the volume (Lebesgue measure) of £1. Eq. (3.2 ) is a Riemann sum generated by partitioning 
into N cells of volume |0|/iV centered at Xj. This yields 

(3-3) ^ = P L 9 ( * ( *' ^ X))A{t ' X)) ^ ,,{X ~ ^ X))dX - 

up to discretization error. Now suppose that the map q(-, X) is invertible for each t, that is X — q~ 1 (t, q). 
Changing the variables in the integral y — q(t, X) we obtain a generic integral approximation 



(3.4) 



9 n = j^| J g{v{t,y),y)Tp ri {x - y)J{t,y) dy, 



where 

(3.5) J= IdetVq" 1 !, 
up to discretization error. 

3.3. Regularized deconvolutions. Define an operator R v by 

Ry[f}(*) = J M*~y)f(y)dy- 

To simplify exposition, suppose that R v is injective. For example, a Gaussian ip^ produces an injective 
operator, which is not difficult to check using Fourier transform and uniqueness of analytic continuation. If R v 
is injective, then there exists the single- valued inverse operator R~ , that we call the deconvolution operator. 
Unfortunately, this operator is unbounded, since R v is compact in L 2 (f2). This is the underlying reason 
for the popular belief that averaging destroys the high-frequency information contained in the microscopic 
quantities. In fact, this information is still there (the inverse operator exists), but it is difficult to recover 
in a stable manner, because of unboundedness. This does not make the situation hopeless, as has been 
recognized for some time. Reconstructing / from the knowledge of R v [f]) is a classical example of an 
unstable ill-posed problem (small perturbations of the right hand side may produce large perturbations of 
the solution) . The exact nature of ill-posedness and methods of regularizing the problem are well investigated 
both analytically and numerically (see, e. g. O [HI HH QTJ [12] ) . Accordingly, we interpret notation i?" 1 
as a suitable regularized approximation of the exact operator. Many regularizing techniques are currently 
available: Tikhonov regularization, iterative methods, reproducing kernel methods, the maximum entropy 
method, the dynamical system approach and others. It is very fortunate that this vast array of knowledge can 
be used for the ODE model reduction. On the conceptual level, our approach makes it clear that instability 
associated with ill-posedness is a fundamental difficulty in the process of closing the continuum mechanical 
equations. 

A family of Landweber iterative deconvolution methods [5], |10j seems to be particularly convenient in 
the present context. In the simplest version, approximations g n to the solution of the operator equation 

(3.6) R V [9]=9 V 
are generated by the formula 

n 

(3.7) g n = Yj!-IQrr, g = 9 r '- 

k=0 

The number n of iterations plays the role of regularization parameter. In (3.7), / denotes the identity 
operator. 

The first three low-order approximations are 

(3.8) g Q = gn n = 0, 

(3.9) 9i=9 n + (I-R v )\9 n ] n=l, 

(3.10) g 2 =^ + (I-R n )m + (I-R v ) 2 m n = 2. 

4. MlCROSCALE PARTICLE CHAIN EQUATIONS 

In this section, the general method outlined above is detailed in the case of one-dimensional Hamiltonian 
chain of oscillators that consists of N identical particles with nearest neighbor interaction. The domain O is 
an interval (0, L). Particle positions, denoted by qj = qj(t), j = 1, . . . , N, satisfy 

< qi < q 2 < ■ ■ ■ < qN < L 

at all times, i.e. the particles cannot occupy the same position or jump over each other. Next, define a small 
parameter 

1 

and microscale step size 

. . L 



The interparticle forces 

(4.2) / ifc =^V(^^ 

are defined by a finite range potential U. We suppose that U'(£) > for all £ within the range. Note that k 
in (4.2 1 can take only two values: j — 1 or j + 1. Also, observe = —fkj, as it should be by the third law 
of Newton, and also that the sign of fjk is the same as sign of qj — q^. This means that the force exerted on 
Pj by say, Pj+x is repulsive. The total interaction force acting on the particle Pj is 

for j = 2,3,...,JV-l. 

Each particle has mass m = M/N — Me, where M is the total mass of the system. Particles have 
velocities denoted by Vj, j = 1, . . . , N. Writing the second Newton's law as a system of first order equations 
yields the scaled microscale ODE system 

(4.3) q 3 = Vj , eMv 3 =f j , j = 1,...,N 
subject to the initial conditions 

(4.4) q ] (0) = q° J , v J (0)=vl 

5. Integral approximation of stresses for particle chains. Mesoscopic continuum 

equations 



In the one-dimensional case stress is a scalar quantity, and (2.101, (2.11) reduce to, respectively, 

N M 

(5.1) a) = -^2-(v j -v^t,x)) 2 ^ v (x-q j ), 



N 

3=1 



and 

N-l 



(5-2) T (Ltf (*' x ) = Y] /;•;• \ {l h- i - Qj) / *Pri( x - s 9j+i - (1 - s)qj)ds. 

j= i Jo 

The sum in (5.2 1 is simplified compared to the general expression, since we have exactly N — 1 interacting 
pairs of particles. 

To obtain integral approximations of stresses, we define interpolants q, v, as in Sect. |3.2| Assuming as 
before that q is invertible and repeating the calculations we get 



L 



M 

(5.3) Tfa(t,x) = -— J (v(t,y)-lf(t,x))^ n (x-y)J(t,y)dy. 

Remark. Many equalities in the paper, including ( |5.3[ ) hold up to a discretization error. To simplify presen- 
tation, we do not mention this in the sequel when discrete sums are approximated by integrals. 
The interaction stress can be rewritten as 

Next we approximate (qj + i — qj)/h by q'(t,Xj). This approximation is in fact exact, provided the interpolant 
is chosen to be piecewise linear. Note also that 

. 1 - 1 

q {t ' X) - " J(t,q(t,X)Y 



Inserting this into (5.4), replacing Riemann sum with an integral and changing variable of integration as in 
Sect. |3.2| we obtain the integral approximation of the interaction stress: 

(5.5) T U <^-/M«)P'H-«>»* 



Equations (5.3 1, (5.5) contain two microscale quantities: J and v. Approximating sums in the definitions of 
the primary averages (2.6), (2.7) by integrals we see that Jf 1 and v* 1 are obtained by applying the convolution 
operator R v to, respectively J and Jv: 



(5.6) 



M 

TP = —R n [J], 



M 

p*T?> = —R v [Jv] 



The discretization error in (5.6) can be made small by imposing suitable requirements on the microscopic 
interpolants. Fortunately, the theory of ill-posed problems allows for errors in the right hand side of integral 
equations. The size of the error determines the choice of regularization parameter. In the present case, 
the error determines the number of iterations needed for the optimal reconstruction, according to so-called 
stopping criteria. These criteria are available in the literature on ill-posed porblems (see e.g. j^]). Detailed 
investigation of these questions is left to future work. 

Denote by TC^\ the iterative Landweber regularizing operators 



R n \ 



k=Q 



(I-Rn) 1 



Applying R^ \ in (5.6) yields a sequence of approximations 



(5.7) 



and a corresponding sequence of closed form mesoscopic continuum equations (written here for an isolated 
system with zero exterior forces) 

(5.8) dtjf>+d x (j?iif) = 0, 



(5.9) 

where Tj 1 , 
(5.10) 



rpf] 

i (int), 



are given by 



■TP- « 

(int), 



o. 



M 

T\ = 

{c),n £ 



(v n (t,y) - v n {t, x)) 2 ip v (x - y)J n (t,y)dy, 



(5.11) 



{int) ,n 



u' 



L 



Jn(t,y)J Jo 



ipr,[x-y 



Jn(t, y) 



ds dy, 



with J n ,v n given by (5.7 1. 



6. Zero-order closure for particle chains 
Let us consider zero-order approximations in detail. The mesoscopic mesh consists of points 

r 



Xp 



1,2 B, 



where B = 1/rj, presumed to be an integer satisfying B -C N. Meso-cells are intervals Ip of length 

L J; = L/B = r)L 1 

centered at xp. 

Suppose that the only primary variables of interest are density p 11 and linear momentum Jpv^ '. These 
variables will be computed by the mesoscale solver. For simplicity, suppose that the meso-solver is explicit in 
time. Then the average density and average velocity will be available at the previous moment of time. Our 
task is to design an update step for computing density and velocity at the next time moment. To construct a 
closed form update step, we need to approximate stress T n in (2.9) in terms o{p r, ,~p n v' 1 . From the knowledge 
of p^iP^v 11 we can approximately recover J and Jv. The zero-order approximation (3.8) corresponds to 

L 



(6.1) 
(6.2) 



J(t,x) 
J(t, x)v(t, x) 



M 
L 

M 



P v (t,x), 
p r 'v 7 '(t,x). 



In other words, the microscale quantities are approximated by their averages. The corresponding closed form 



approximations for stress are obtained by inserting (6.1|, (6.2) into (5.10), (5.11): 



(6.3) 



Q (t,x) 



(0"(t, y) - v\t, x)Y %{x - y)]F{t, y)dy, 



(6.4) 



(mt),0 



U' 



M 



shM 



V 



L]P{t,y) 



ds dy. 



For computation, a numerical quadrature should be used. In this regard, note that all average quantities are 
computed on the mesoscale mesh, while the formulas (2.10), (2.11) are fine scale discretizations. Therefore, 



one might wonder if a straightforward mesoscale quadrature of (6.3), (6.4) is too crude. A better approach 
is to interpolate p v , I? 1 ' by prescribing approximate particle positions q~j and velocities ijj, j = 1,2, ...,N, 



compatible with the given p ri ,v ri . Once this is done, (6.3), (6.4) can be discretized on a fine scale mesh with 
mesh nodes cjj . 

Interpolants cannot be unique. For zero-order closure, we are choosing positions and velocities that pro- 
duce the given average density and average velocity. Clearly, there are many different position-velocity 
configurations with the same averages. The choice made in the paper is motivated by the practical require- 
ment of achieving low operation count, as well as by certain expectations about the nature of dynamics. 
From the continuum mechanical point of view, if a system can be adequately modeled by balance equations 
of mass and momentum, then it must have have a trivial energy balance. Most often this means that the 
deformation is nearly isothermal. To mimic such isothermal dynamics we suppose that at each time step, 
there exists a positive number n 2 (it can be called "upscaling temperature") such that 



(6.5) 



(uj - ytfax^ftp^xp - qj) = k 2 



Here the summation is over all particles located in a meso-cell I p. The temperature k 2 is the same for all 
= 1,2,..., B. We emphasize that the actual value of k 2 is not as important as the fact that its value 
is the same for all meso-cells. This is because (6.5) would yield constant mesoscale mesh node values of 
Q As a result, the finite difference approximation of d x T? e s on the mesoscale mesh is identically zero. 
We interpret this by saying that convective stress does not contribute to the mesoscopic dynamics in the 
isothermal case. Another observation is that k need not be the same at different moments of time, so our 
assumption is somewhat more flexible than the standard isothermal deformation approximation. Also, we 



note that validity (6.5) depends on the choice of i]. For bigger 77, it is more likely that (6.5) holds for the 
same underlying microscopic dynamics. Details on this are provided below in Section |6.2| Additionally, 
other features of the microscopic dynamics should be taken into account. Most importantly, interpolated 
■ v v (qj) must be such that the collection qj,ij,j — 1, 2, . . . , N conserves microscopic energy £: 



velocities 



(6.6) 



1M A 2 

3=1 



■U{Q) 



where U(Q) is the microscale potential energy corresponding to the positions qj. 

6.1. Prescribing particle positions. The objective of this section is to assign approximate particle posi- 
tions qj. We start by interpolating J. The simplest interpolant is piecewise-constant: J(t, x) ps Ylf=x ^(*; x p)xp{ x ) 

Ss=i ~M~P v (ti x fi)Xl3{ x ) where x/3 is the characteristic function of the meso-cell Ip. A simple choice of the 
position map compatible with this interpolant is a piecewise linear map having the prescribed constant value 
of J in each meso-cell. In practical terms, this means that in each meso-cell, particles are spaced at equal 
intervals from each other. The local interparticle spacing 



(6.7) 



M 



is determined by the mesh value of the average density. To explain (6.7), note that the total mass of particles 



contained in the meso-cell Ip can be approximated by p n {t 1 xp)L^. Dividing by the mass M/N of one particle, 
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we obtain an approximate number of particles inside I, 



np = p v (t,xp)I v 



N 



and thus A.p = L v /np. We emphasize that qj are chosen based only on the known mesh values of the density 
p v , and that qj will be different from the actual particle positions qj. 

Now we approximate the integral in (6.4) by its Riemann sum generated by the partition {q~j,j = 
1,2,. ..,N}: 

N-l 



(6.8) 



(int),0 



JV-1 .1 

J2 U' (N(q j+1 - qjj) (q j+1 - qj) \ 

.7 = 1 J ° 



ip v (x - sqj+i - (1 - s)qj) ds. 



6.2. Prescribing particle velocities. In order to approximate the convective stress in the fine scale dis- 
cretization of (6.3 1, we need to choose approximations Vj of the true particle velocities Vj. The choice of Vj 
must satisfy (6.6) and be compatible with the available average velocity at the mesoscale mesh nodes. 
For each qj € Ip, we set 

Vj =v v fj + 5v?, 

where is the local average velocity, and Svj is a perturbation to be defined. 

Next, we show that the energy-conserving collection of Svj velocity always exists, provided its upscaling 
temperature is suitably prescribed. This prescription will be based only on the available mesoscale informa- 
tion. For definitiveness, in the rest of this section we suppose that ip v satisfies the following condition: 

(6.9) ip v (xp - qj) > 0, if qj G Ip. 

To make the algebra simpler, we make another assumption: for each (3 — 1,2, ... ,B, 

N 

(6.10) f(vj)*Pri(xp ~ qj) « /(%)^( x /3 - 



j=i 



where / is either Vj or (vj) 2 - The second summation is over all j such that qj £ Ip. Assumption (6.10) holds 
when iprj(xp — y) is small outside of I p. 

Averaging of Vj should produce the known average velocity til This yields an equation for 8vj\ 



(6.11) 
Since 



M 
iV 



Y Vjtp n {xp - qj) 



J p- 

—7] —71 

Pp v p- 



M 
~N 



2 Vjipni 



xp - q 3/ 



jeJp 



53 ^( X P - Qj) + ^7 Yl 5v fav( x P - Qj) 



M 
~N 



N 



(6.11) holds provided 



M 

iV 



(6.12) 

Now we look for perturbations in the form 
(6.13) 



53 8v?ip v (i 



'/9 - Qj) 



5v? = 



^{xp-qj)'' 



where the are to be determined. Next, narrow down the choice of a? by setting 



(6.14) 



where dj is either one or negative one. To satisfy (6.12) we need rip to be even (one more point can be easily 
inserted if the actual np is odd); in addition, the number of positive and negative d\j must be the same. 
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To simplify further calculations, wc write conservation of energy (6.6) in the form 
(6.15) 



where Kg are any numbers satisfying 



(6.16) 



Kg > 0, K f> 



2N 



8=1 



Our goal now is to show that there is a choice of Kg, k and t such that Sv^ defined by (6.131, (6.14) 
satisfy equations (6.5), (6.15), and ( 6.16[ ). Inserting (6.13) into (6.5) and (6.15) yields, respectively, 

1 



(6.17) 
(6.18) 



* 2 E 



^r,( X 8 - Qj) 



< 2 E 



1 



{■*l> v {x B -Qj)) 2 



Combining these equations we get 
(6.19) ? = ^ ( £ 



^{xg - qj) 



(6.20) 



5.2 1 



Substituting into (6.16) yields the choice of k: 



E 



1pr,( x 9 - Qj) 



= Ka 



(6.21) 




£-U{Q) 



2 JV ^ 

(9=1 



(^) 2 «/3 



The choice of all constants now should be made as follows: 
1) Given £,lfl,qj, find k by (6.21); 



2) Determine Kg from (6.20); 



3) Determine t from (6T9 ) 



4) Choose 5vf by (16.13 



6.14) 



3=1 t-'jZJp il>r,(.xp-qj) 



Note that step 4 introduces non-uniqueness, but we are concerned only with existence of suitable velocity 
perturbations. The actual choice of Sv^ will not change the mesoscopic discretization of momentum balance 



equation. Indeed, once qj,Vj are chosen, we can approximate the integral in (6.3) (for x at the mesoscale 
mesh nodes) by a Riemann sum corresponding to the partition qj of (0,L): 



(6.22) 



-^(c),o(*> Xa ) ~ 



b=i je J,) p 

~ E (tiv?) 2 i/J v (x a - Qj) 
je,J a 



a = 1,2, . . . ,B. 



Qj)rf 



Therefore, the mesoscale mesh values of T?* Q are all equal. This implies that a finite difference approximation 
of d x T^ Q on the mesoscale mesh vanishes. The conclusion is that for isothermal dynamics, any suitable 
choice of a velocity perturbation produces a convective stress that has zero divergence on the mesoscale. 

li 



6.3. Zero-order isothermal continuum model. Combining the approximation d x T 7 ^ Q = with (|6.3|), 



(|6.4|) we obtain an isothermal zero-order continuum model 

dtp 11 + 0x(p"w") 



(6.23) 
(6.24) 



= 0, 
= 0, 



where T^„ i ^ n is given by an integral expression (6.4 1 (or by a discretization (6.8|). We can interpret 



- (int),0 



interaction stress as pressure. Then (6.4) provides dependence of pressure on density. 



which is non-local in 



space and non-linear. For small h (large N), it can be approximated by 



M 



ipr, (x - y) dy, 



jo U> \p n (t,y 

which is still non-local. In the limiting case 77 — >• 0, observing that convolution with is an approximate 
identity, this equation can be reduced a local equation of state 



(int),0 



-u' 



M 



P v (t,x) 



If p n is nearly constant, this equation can be linearized to produce a classical gas dynamics linear equation 



of state. This shows that zero-order closure (6.4) generalizes several classical phenomenological equations of 



state. The connection between micro- and mesoscales is made explicit in (6.4). Using higher order closure 



approximations, one can obtain other non-classical continuum models worth further investigation. 



7. Computational results 

In this section, the method developed in the previous sections is tested for a chain of N = 10, 000 to 
N = 80, 000 particles interacting with a non-linear finite rage potential 

0, 



(7.1) 



U(0 = 



2-p 



if £ e (0,x*] 
if £ > x+ 



where p > 1, x* = aL, awl and C r is material stiffness. This potential mimics a Hertz potential used 
in modeling of granular media. Particles in this model are centers of lightly touching spherical granules 
arranged in a chain, and the ODEs model acoustic wave propagation in this chain. The microscale equations 



are (4.3), (4.4) with initial conditions given below. The forces include the pair interaction forces defined by 



U, and the exterior confining forces acting on the first and last particles. The parameters of U are chosen so 
that all particles stay within the interval [0, L] for the duration of a simulation. To ensure that particles do 
not leave the interval [0, L], its endpoints are modeled as stationary particles that interact with the moving 
particles with forces generated by the same potential U. If needed, stiffness of walls can be increased by 
using a different value of C r . 

Let xp he the centers of the mesocells Ip, (3 = 1,2,..., B, as defined in Section [6] Next, let a window 
function ^p(x) be the characteristic function of the interval [—^r/L, \rjL). The average density and momentum 
are defined by (2.6) and (2.7), respectively, and the average velocity is 



Efci^- 0M£/3 - q,{t)) 



The average density p v evaluated at the center xp of a, mesocell Ip can be written as 



(7.2) 



p}(t)=pT>(t,Xp) 



3=1 



V 



N L 



that shows that average density is proportional to the scale separation B/N. 
We solve microscopic equations (4.3), (4.4) subject to the initial positions 



<l;, = \ 3 ~ h, j = 1, 2, . . . ,N, h=jj, 
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FIGURE 1. N = 40,000, B — 50. Dashed line: Jacobian J(£,xg); solid line: its mesoscale 



approximation jjp v (t, xp), ft — 1,2, . . . , B according to (6.1 1. Blowup of results at t = 0.01 
shows discrepancy between J(t,xp) and jj~p n {t,xp). 





t=0.03 



0.5 




exact interaction stress T^ int ^(t, xp); solid 



FIGURE 2. N = 40,000, B = 50. Dashed line: 
line: its approximation ^ (i, xp), ft = 1,2, ... ,B, defined in ( 6.4 ) . Blowup of results at 
t = 0.001 shows difference between exact stress and its approximation. 



and the initial velocities 



7: 

7(" 
0, 



5 



if < 9° < § , 
if f <<Z°<f, 
if f < q° < L 



using the Velocity Stormer-Verlet method. We use L = 1, p = 2, a = 1, 7 = 0.3 and C r = 100. This velocity 
profile initiates an acoustic wave that propagates to the right. The stiffness constant C r can be used to 
ensure that particles have only small displacements from their equilibrium positions. Using initial velocity 
with higher 7 would require a higher value of C r to enforce smallness of typical particle displacements. 
With fixed N = 40,000, B = 50, we integrate microscopic equations (4.3), (4.4) until the acoustic wave 



reaches the right wall, interacts with it and is about of being reflected to the left, which corresponds to 
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t = 0.07. To capture the most interesting dynamics, we present snapshots of results at times t — 0, 0.01, 
0.03, 0.05, 0.06 and 0.07. To test our closure method, we compute microscopic positions qj and velocities 



Uj, 



i = i, 



, N, at every time step and use them to evaluate primary mesoscopic variables: average density 



Pp and average velocity x fL, a t mesocell centers xp, f3 = 1, . . . , B. These mesoscopic quantities are defined 
in (2.6), (2.7 1 (see also (|5.6| ). They are then employed in computing of the zero-order approximation 
T^ n( j (t, xp) defined in (|6.4|). We compare this mesoscopic approximation with the "exact" microscopic 
interaction stress Tjl_^ (t,xp) defined in (5.2), and also test other approximations given by (6.1), (6.2). 



' (int) ' 

Comparing Vj, j = 1, 



,N and v}, (3 



1, 2, . . . , B (not shown here) we find that micro- and mesoscale 



velocities are essentially indistinguishable during the simulation time. 
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FIGURE 3. N = 40,000, B = 50. Convective stress T?Jt,xp), /? = 1,2,..., B, defined in (6.3) 



In Fig. [Tj we analyze microscopic Jacobian J(t, Xp) together with its zero-order mesoscopic approximation 



jT]f(t,xp) obtained according to (6.1). In this and other figures, we plot "exact" microscopic quantities 
using a dashed line while mesoscopic quantities are depicted with a solid line. Results shown in Fig. [I] 
indicate that jjp n (t,xp) exhibits some oscillations whose amplitude is about 10~ 3 as compared to Jacobian 
J(t, xp). The oscillations are likely caused by the choice of a window function tp. For computational testing, 
we chose ip to be a characteristic function. The main reason was to try "the worst case scenario" concerning 
smoothness of ip. We expected that this window function would produce more oscillations than a smoother 
ip. A good agreement between our approximation and the direct simulation results strongly suggest that 
the proposed method is viable. We believe that it should perform better with a smoother choice of ip. The 
oscillations present in j^~p r ' '(t , x p) are amplified in the approximated stress T? int \ o(f,^), due to the rather 
high stiffness constant C r = 100, as shown in Fig. [2j We al so co mpare microscopic J(t,xp)v(t,xp) with 
its zero-order approximation jjJP^jX^v^ftjXp) according to (6.2). Graphs are not shown here but we find 



that these quantities agree very well similar to micro- and mesoscale velocities. This is expected since the 
average density is approximately identity with small oscillations. Finally, we verify that the dynamics is 
quasi-isothermal by plotting the convective stress T^(t,xp) defined in (5.11) in Fig. |3j As can be seen, 
fluctuations in the convective stress do not exceed 10 -4 throughout computational time, therefore, the kinetic 
energy of velocity fluctuations is small. 

We next tested the effect of the scale separation on the quality of the zero-order approximation. With 
fixed B = 50, we allowed iV vary from 10, 000 to 80, 000 and followed the evolution of mesoscale quantities of 
interest: ■^p n (t, xp) and TJ^^ (t, xp). Snapshots of these functions at the same representative time t — 0.01 
are plotted in Figs. [|]and[5j respectively, with N = 10, 000, N = 20, 000 and N = 80, 000. The results with 
N = 40, 000 at the same time are given the middle top panels in Figs. [T] [2] for comparison. It is clear that as 
scale separation increases, oscillations in both ^p v (t, xp) and 1^ int \ (i> x p) diminish and when N — 80, 000, 
the exact microscopic quantities and their mesoscale approximations are almost indistinguishable. 
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Figure 4. Effect of the scale separation on jjp n - B = 50 is fixed, N varies, data is taken at 
the same t = 0.01. Dashed line: Jacobian J(t,xp)\ solid line: its mesoscale approximation 
±7f>(t,x f3 ),p = l,2,...,B. 
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Figure 5. Effect of the scale separation on TV. 
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B = 50 is fixed, N varies, data is 



taken at the same t = 0.01. Dashed line: exact interaction stress T? int \ (t, xp); solid line: its 
mesoscale approximation T? int -. Q (t, xp), j3 = 1, 2, . . . , B. 
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Figure 6. Example with imposed high frequency oscillations. N — 10,000, B = 50. 
Dashed line: exact velocity Vj, j — 1,2, ...,7V; solid line: average velocity vp, (3 — 
B. 



1,2, 



2r 

1.5 
1 

0.5 
0. 



t=0 



2 
1.5 



0.2 



0.4 



0.6 



0.8 











t=0.01 



Figure 7. Example with imposed high frequency oscillations. N = 10, 000, B = 50. 
Dashed line: Jacobian J(t,xp); solid line: its mesoscale approximation j^fp^^xp), j3 = 
1, 2, . . . , B (compare with Fig. [T]). 



In the above example, fluctuations of microscopic velocities about their average values were very small and 
the zero-order approximation worked well. Next we show that if microscopic velocities have high fluctuations 
then the zero-order approximation is not capable of captioning an appropriate dynamics. 
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Wc demonstrate this by imposing high frequency k oscillations with relatively large amplitude a on the 
nonzero portion of the initial velocity used in the previous experiments. The initial velocity is 

f J + asin(^q°), if 0<g?< f, 

v°=l 7 (-H° + 2)+asin(^), if § < q o < f , 

1 0, if f < <?° < L. 

and it is plotted in the left panel of Fig. [6] We use a — 5 and fc = 20 that gives one period of imposed 
oscillations per mesocell. This microscopic initial velocity has a property that the average velocity at time 
t = is exactly the same as in the previous example. Simulations were done with N = 10, 000 until the 
same t — 0.07. 




Figure 8. Example with imposed high frequency oscillations. N = 10,000, B = 50. 
Dashed line: exact interaction stress T^ int ^(t, xp); solid line: its mesoscale approximation 
T (int),o& x P) (compare with Fig. [2j). 

The right panel of Fig. [6] shows a typical microscopic velocity profile together with its average velocity 
(taken at t = 0.01): to the left from the wave front, the microscopic velocity has large frequency oscillations 
(due to dispersion?) with an amplitude sometimes exceeding the initial amplitude by a factor of 1.5 and to 
the right from the wave front, the microscopic velocity is zero. Clearly, the average velocity is very different 
from the microscopic velocity. 

Analysis of microscopic Jacobian J(t,xp) and mesoscopic jr]f(t,xp) reveals that these functions have 
qualitatively the same dynamics as micro- and mesoscale velocities, respectively, shown in Fig. [6] We plot 
the former in Fig. [7] where the left panel has graphs at t = while the right panel shows typical structure 
with data taken at t = 0.01. It is interesting to note that the zero-order approximation 1? int \ q{PiX^) to 
the interaction stress T^ int ^(t,xp) plotted in Fig. pmoes not agree in those areas that were affected by large 
magnitude oscillations in microscopic velocities while agrees well in those areas to which oscillations have 
not come yet. This finding suggests that indeed the zero-order approximation should not be used for large 
frequency oscillations in microscopic velocities and a higher order approximation is needed. 

Finally, in Fig. |9| we plot the convective stress T^(t,xp) whose large values confirm that oscillations 
in microscopic velocities are much bigger during the computational time than those in the first example. 
When the initial velocity has fluctuations with frequency higher than k — 20, discrepancy between micro- 
and mesoscale quantities is even more pronounced. 

8. Zero-order closure: applicability 

Zero-order closure is very similar to the use of the Cauchy-Born rule in quasi-continuum simulations of 
solids. Here, the nodes of the mesoscale mesh can be thought of as "representative particles" . These particles 
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Figure 9. Example with imposed high frequency oscillations. N = 10,000, B = 50. 
Convective stress T^(t, xp), j3 = 1, 2, . . . , B (compare with Fig. |3j). 



are moved with the average velocity, while the velocities of other particles are assigned by interpolation. A 
construction of an intcrpolant should take into account the physics of the microscopic model such as energy 
conservation. In the computational example of Section [7J zero-order approximation turns out to be quite 
accurate, when non-oscillatory initial conditions are imposed. In this case, we found that approximate and 
exact stresses agree rather well, and this agreement becomes better with increasing scale separation. 

This does not mean that zero-order closure always works well. Our numerical simulations suggest that 
applicability of zero-order closure is determined by initial conditions, exterior forces and interaction potential 
(arranged in order or importance). 

Approximating functions by their averages we neglect fluctuations. Therefore, the initial velocities should 
have small fluctuations. Initial positions should be chosen so that the number of particles in a meso-cell 
varies slightly from one cell to another. The initial velocity fluctuations in our first example are small, and 
convective stress at later times is by three orders of magnitude smaller than interaction stress. This remains 
true on the time interval sufficient for the traveling wave to reach the opposite end of the chain. 

For one-dimensional problems, convective stress is proportional to the kinetic energy of velocity fluctua- 
tions. This kinetic energy can be naturally associated with upscaling temperature. Relative smallness of the 
convective stress mens that upscaling temperature is nearly zero. Therefore, the corresponding dynamics 
can be termed cold. We also note that cold dynamics is a special case of isothermal dynamics, considered 
in Section [6] As has been remarked earlier, isothermal dynamics implies that divergence of the convective 
stress is nearly zero on mesoscale, and thus can be neglected compared with the divergence of the interaction 
stress. 

Another consideration is related to inhomogeneity in actual particle distribution. In our example, devi- 
ations of about 4% in relative particle positions produced visible oscillations in the approximation of the 
interaction stress. This amplification of small perturbations is due to the stiffness of the interaction poten- 
tial. However, the same stiffness prevents particle aggregation, keeping the interparticle distances bounded 
from below. Bounds from above are difficult to enforce with the chosen potential because it does not have 
a potential well. The isolated particle system with this potential would just fall apart. This phenomenon is 
common place for granular materials. The particles remain confined to the domain (container) only because 
they are repelled by the walls. Walls have very little direct influence on the interparticle distances in the 
systems' interior. Therefore, applicability of zero-order closure also depends on the stiffness of the problem, 
and more generally on how well the potential enforces uniform particle distribution. In that sense, zero- 
order closure makes a reasonable approximation for lattice systems modeling small deformation of solids at 
constant temperature. 
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To further understand limitations of zero-order closure, consider the effect of increasing the order n of the 
Landwcbcr approximations (3.7 1. The Fourier transform the kernel of I — R v is equal to 



1 - er v n 

It is very small for £ close to zero, and then increases to one as |£| goes to infinity. Therefore, I — R v acts as 
a filter damping low frequencies and thus emphasizing higher frequency content of the signal. Higher order 
approximations amount to applying convolutions Ysk=x(I~ ^v) k t° mesoscale averages. As n increases, high 
frequency content of the reconstruction will be increasingly amplified. This suggests that systems capable 
of producing large fluctuations should be handled with higher order approximations. 

A related comment is that averages of fluctuations can become additional state variables in a mesoscale 
continuum model. A familiar example is the use of the averaged energy balance equation (see [T] for deriva- 
tion), in addition to the mass and momentum balance. The energy balance equation describes evolution of 
the density of kinetic energy of velocity fluctuations. An intriguing question here is how the model with just 
two equations of balance but high order closure approximation compares with a zero-order closure model 
containing all three balance equations. In classical physics, additional balance equations are often introduced 
as a means of compensating for errors introduced by replacing state variables with their averages. Use of 
higher order closure could offer an alternative to this approach. Indeed, suppose that one is interested only 
in tracking density and velocity on mesoscale. The corresponding two balance equations contain only two 
microscale quantities: velocity field v and the Jacobian J of the inverse position map q . If v and J can be 
accurately reconstructed from their averages, we do not need to deal with the energy balance equation. This 
observation offers a new way of reducing computational cost. Higher order approximation are more expen- 
sive than zero-order, but using more balance equations also increases computational cost. We also note that 
increasing the order of closure approximations involves repeated convolutions with the window function ip^. 
On the other hand, simulating an energy balance involves numerical integration of an additional non-linear 
integral-differential equation, a much more difficult task. 



9. Conclusions 

We propose a closure method that gives closed form approximations for mesoscale continuum mechanical 
fluxes (such as stress) in terms of primary mesoscopic variables (such as average density and velocity) . Our 
closure construction is based on iterative regularization methods for solving first kind integral equations. 
Such integral equations are relevant because mesoscopic density and velocity are related to the corresponding 
microscopic quantities via a linear convolution operator. The problem of inverting convolution operators is 
unstable (ill-posed) and requires regularization. Use of the well known Landweber iterative regularization 
yields successive approximations, of orders zero, one, two and so forth, to interpolants of particle positions 
and velocities in terms of available averages. Closure is achieved by inserting any of these approximations 
into the equations for fluxes instead of the actual particle positions and velocities. Low order approximations 
are simpler to implement, while higher order approximations can be used to more accurately reproduce the 
high frequency content of the microscopic quantities. 

The above general strategy is applied in the paper to spatially averaged dynamics of classical particle 
chains. We focus on the simplest zero-order approximation and show numerically that it works reasonably 
well as long as initial conditions have small velocity fluctuations. The case of large fluctuations in velocities 
should be handled by higher order approximations. 
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